
from rasterutils import checkRasterSupport, layerAsArray, writeGeoTiff

import gdal

# let's try it with numpy:
import numpy
import math

# this lets us use function-like operators:
import operator

def I(x):
    """an identity operator"""
    return x

# a list of (name, function) tuples.
# functions must take one argument only.
# sorry, but them's the rules
# its an array not a dict, because I want to control the order
# they appear in the dropdown.

_funcs = [
    (' ', I),
    ('sqrt', numpy.sqrt),
    ('sin', numpy.sin),
    ('cos', numpy.cos),
    ('tan', numpy.tan),
    ('exp', numpy.exp),
    ]

# a list of (name, operator) tuples.
# 
_ops = [
    ('+',operator.add),
    ('-',operator.sub),
    ('*',operator.mul),
    ('/',operator.div),
    ('==',operator.eq),
    ('>',operator.gt),
    ('<',operator.lt),
    ('>=',operator.ge),
    ('<=',operator.le),
    ]


def getFunc(name):
    """ search through the functions to find one with a name"""
    for f in _funcs:
        if f[0] == name:
            return f[1]
    print "Function not found... returning identity"
    return I


def getOp(name):
    """ search through the operators to find one with a name """
    for o in _ops:
        if o[0] == name:
            return o[1]
    print "Operator not found... returning 'plus'"
    return operator.add

# return list of operator names
def ops():
    return [x[0] for x in _ops]


def funcs():
    """ return list of function names """
    return [x[0] for x in _funcs]


def createRaster(filePath, formula):
    """ create a raster from a formula and stick it in a file"""

    # must have one layer being used (otherwise we've no extent)

    if not formula['useLayer1'] and not formula['useLayer2']:
        print "Must have at least one layer selected..."
        return

    # compute the extent of the result, and work out the
    # subjects of the left and right functions
    
    # if using both, must be roughly same extent
    if formula['useLayer1'] and formula['useLayer2']:
        extent = checkRasterSupport(formula['layer1'],formula['layer2'])
        if not extent[0]:
            print "Layers not close enough in grid size/posn"
            return
        extent = extent[1]

        leftSubject = layerAsArray(formula['layer1'])
        rightSubject = layerAsArray(formula['layer2'])
        
    else:
        if formula['useLayer1']:
            layer = formula['layer1']
            leftSubject = layerAsArray(formula['layer1'])
            rightSubject = float(formula['constant2'])
        else:
            layer = formula['layer2']
            leftSubject = float(formula['constant1'])
            rightSubject = layerAsArray(formula['layer2'])
            
        le = layer.extent()
        extent = [le.xMin(),le.yMin(),le.xMax(),le.yMax()]

    # then its as simple as this:
    result = getOp(formula['op'])(
        getFunc(formula['func1'])(leftSubject),
        getFunc(formula['func2'])(rightSubject)
        )

    # which is f1(leftSubject) %op% f2(rightSubject)

    # now make a geotiff out of it.

    
    writeGeoTiff(result, extent, str(filePath))

    # and we're done
    return

    
    
